Unveiling interactions between intervertebral disc morphologies and mechanical behavior through personalized finite element modeling

Introduction: Intervertebral Disc (IVD) Degeneration (IDD) is a significant health concern, potentially influenced by mechanotransduction. However, the relationship between the IVD phenotypes and mechanical behavior has not been thoroughly explored in local morphologies where IDD originates. This work unveils the interplays among morphological and mechanical features potentially relevant to IDD through Abaqus UMAT simulations. Methods: A groundbreaking automated method is introduced to transform a calibrated, structured IVD finite element (FE) model into 169 patient-personalized (PP) models through a mesh morphing process. Our approach accurately replicates the real shapes of the patient's Annulus Fibrosus (AF) and Nucleus Pulposus (NP) while maintaining the same topology for all models. Using segmented magnetic resonance images from the former project MySpine, 169 models with structured hexahedral meshes were created employing the Bayesian Coherent Point Drift++ technique, generating a unique cohort of PP FE models under the Disc4All initiative. Machine learning methods, including Linear Regression, Support Vector Regression, and eXtreme Gradient Boosting Regression, were used to explore correlations between IVD morphology and mechanics. Results: We achieved PP models with AF and NP similarity scores of 92.06\% and 92.10\% compared to the segmented images. The models maintained good quality and integrity of the mesh. The cartilage endplate (CEP) shape was represented at the IVD-vertebra interfaces, ensuring personalized meshes. Validation of the constitutive model against literature data showed a minor relative error of 5.20%. Discussion: Analysis revealed the influential impact of local morphologies on indirect mechanotransduction responses, highlighting the roles of heights, sagittal areas, and volumes. While the maximum principal stress was influenced by morphologies such as heights, the disc's ellipticity influenced the minimum principal stress. Results suggest the CEPs are not influenced by their local morphologies but by those of the AF and NP. The generated free-access repository of individual disc characteristics is anticipated to be a valuable resource for the scientific community with a broad application spectrum.


CONSTITUTIVE MODELING OF THE IVD
The IVD material model considers (1) a solid phase comprising structural macromolecules such as collagen, elastin, and proteoglycans, alongside cells and (2) a fluid phase consisting of water and solutes (Malandrino et al., 2015).The biphasic-swelling (BS) theory, as detailed by Mow et al. (1980Mow et al. ( , 1989)), delineates both the equilibrium and transient mechanics of charged soft tissues in IVDs.This theory presents each tissue as a composite material featuring a charged solid porous phase saturated by interstitial fluid, thereby enabling the simulation of fluid pressurization and movement within the disc.
This study characterizes the behavior of the entire disc through an osmo-poro-hyper-viscoelastic model.This comprehensive model integrates the constitutive tissue of the bony endplate, treated as a linear poroelastic material (Malandrino et al., 2011) through shell elements.Furthermore, the annulus, nucleus, cartilage, and transition zone, represented with second-order hexahedral elements, employ the BS model used to simulate poromechanical interactions within a poro-hyperelastic matrix saturated with intra-and extra-fibrillar fluid (Wilson et al., 2005b), including the Donnan osmotic pressure gradient effects (Urban and Maroudas, 1981), which, according to Darcy's law, governs the flow of interstitial fluids through the tissue's porous matrix (Wilson et al., 2005b).In addition, the model considers viscoelastic collagen fibers present only in the AF Wilson et al. (2006a) as rebar elements.
The total stress tensor σ tot is expressed as the superimposition of the effective stress σ eff (defined in 1.0.1) of the solid skeleton within the pores, a fluid pore pressure component p, and Darcy's law: where I is the identity tensor, q is the fluid mass flow to the spatial gradient of pore pressure ∇p, and κ is the hydraulic permeability tensor of the tissue.Also, the fluid flow can be expressed by: where u f is the pore fluid velocity, and n f represents the total water fraction, i.e., the porosity of the medium.
Due to the fixed charges, the cation concentration inside the tissue is higher than in the surrounding body fluid (Wilson et al., 2005b).This excess of ion particles within the matrix creates the Donnan osmotic pressure, ∆π, which drives the fluid flow, causing the swelling of the tissue (Urban et al., 1979).
Incorporating the osmotic pressure into Equation S1, where Schroeder et al. (2007) adapted this equation to the IVD, the hydrostatic fluid pressure p is defined as: u w is the water chemical potential, linked with the pore pressure degree-of-freedom generated by the interstitial fluid permeation effects through the permeability (introduced in 1.0.5) by applying Darcy's law to describe a relationship between fluid flow and the swelling pressure.Therefore, fluid flow between the different tissues of the model depends on the tissue-specific mappings of permeability.∆π represents the osmotic pressure gradient generated by the difference between the internal and external salt concentrations (more details in Section 1.0.3).

Solid matrix -Non-fibrillar Part
While it is a foundational assumption that the organic solid matrix and the interstitial fluids are intrinsically incompressible, the overall solid matrix displays compressibility attributable to its porous structure (Mow et al., 1980).Addressing this, Wilson et al. (2005b) have introduced a formulation specifically designed for large deformations.This is particularly relevant as substantial pore collapses can generate incompatibilities between the hypotheses of solid matter incompressibility and porous solid compressibility.Consequently, they propose that the effective bulk modulus should be considered strain-dependent.To model the behavior of the non-fibrillar solid matrix within this framework, Wilson et al. (2005a) adopted the following compressible neo-Hookean model: where J is the determinant of the deformation gradient tensor F , and I 1 is the first invariant of the left Cauchy-Green strain tensor B = F • F ⊺ .The bulk, K m , and the initial shear modulus, G m , are defined as: where E m is the Young's Modulus and ν m the Poisson's ratio.
For a solid fraction of 1, if the fluid is fully expelled, there are no pores, making the entire matrix incompressible.Conversely, as the solid fraction goes to zero, the volume fraction of the pores goes to 1, rendering the solid matrix fully compressible.To include this in the model, Wilson et al. (2005b,a) proposed the following expression for the Poisson's ratio: where n s,0 and n s are the initial (in the unloaded and non-swollen state) and current solid volume fraction.
Note that ν m goes to 0.5 for a solid fraction of 1 and becomes zeros when the solid fraction goes to zero.This term of Equation ensures that tissue compressibility vanishes as J approaches n s,0 , i.e., the fluid is fully expelled.By substituting Equation S8 into Equation S6 and S7, and by eliminating E m we get the following expression for K m : Because the second law of thermodynamics needs to be fulfilled, the constitutive relation for Km cannot be directly substituted in Equation S5.Therefore, it is first implemented into the energy function and derives the new Cauchy stress from the energy function: where C is the Cauchy-Green deformation tensor.Now, the macroscopic stress-strain response of the solid matrix is determined by G m , n s,0 , and the current deformations of the homogenized poroelastic continuum.This response follows the Cauchy stress of the non-fibrillar matrix to describe the material's finite strain behavior, as Wilson et al. (2005bWilson et al. ( ,a, 2006a,b) ,b) detailed initially and then adapted by Schroeder et al. (2008) for the IVD: 1.0.2Annulus Fibrosus Collagen Fibers -Fibrillar Part In the AF, collagen fibers exhibit a unidirectional viscoelastic mechanical response.This behavior is modeled by incorporating finite strains in a Zener viscoelastic model with two non-linear springs.Assuming that the fibrils only resist tension, the Cauchy fibril stress tensor in a unit area for viscoelastic fibrils (Wilson et al., 2006a,b) can be expressed as: where ψ is the elongation of the fibril, P f is the first Piola-Kirchhoff fibril stress, and ⃗ e f is the current fibril direction.
The mechanical behavior is further refined by integrating an elastic component, represented by an exponential stress-strain relationship, with a viscoelastic component (Wilson et al., 2006a).The collagen fibrils' viscoelasticity is described as a strain-dependent (Charlebois et al., 2004) behavior with the following two-parameter exponential stress-strain relationships of the Zener model: and for ε e ≤ 0 (S14) Here, E 1 , E 2 , k 1 , and k 2 are positive material constants; ε f represents the total fibril strain, while ε e is the strain in the spring s 2 that interacts in series wit the dash-pot, η.The transient stress includes the dashpot effect that partially relaxes P 2 .At equilibrium (i.e., full relaxation), P 2 vanishes completely.
The derivation of P 2 as a function of the fibril strain ε f is given in Wilson et al. (2006a).

Pressure Component -Osmotic swelling
Characterizing the behavior of the fluid by Darcy's law alone is insufficient.Its response depends on tissue deformations and pressure gradients and is also affected by the chemical fixed charge densities that create swelling driving forces.In the BS model, under free swelling conditions, the volumetric expansion results in an increase in the Jacobian determinant.
The Donan osmotic potential describes swelling behavior (Malandrino et al., 2015), assuming that electrolyte flux can be neglected in mechanical studies of charged materials.Accordingly, the internal and external osmotic pressures are represented by the classical Van't Hoff equation (Huyghe and Janssen, 1997), and assuming that the osmotic components are instantaneously equilibrated with the external bath, the osmotic pressure gradient ∆π is given by (Wilson et al., 2005b): where R is the gas constant, and T is the absolute temperature.The internal and external osmotic coefficients ϕ int and ϕ ext multiply the terms related to concentrations of mobile cations and anions, respectively.The average of the internal and external activity coefficients of the ions is represented by γ ± int and γ ± ext (assuming γ ± = γ + γ − ).These osmotic and activity coefficients were implemented as Huyghe and Janssen (1997); Huyghe et al. (2003) proposed and as other authors such as (Wilson et al., 2005b), (Schroeder et al., 2007) and Galbusera et al. (2011) adopted in their respective studies.The external concentration of salt and the proteoglycan fixed charge density are denoted by c ext and c f,exf , respectively.

Tissue model parameters and relation to composition measurements
To elucidate the relationship between biphasic/poroelastic models and IVD deformations quantified by J, it's essential to connect the variables in the equations for the non-fibrillar solid matrix's effective stress S11 and the osmotic potential S15.The proteoglycan fixed charge density is determined by the ratio of the normal fixed charge density (c f ) in milliequivalents per milliliter of total fluid to the extra-fibrillar water (n f,exf ), as defined by Ruiz et al. (2016): where n f,exf is derived as: In this Equation, φ ci indicates the intrafibrillar water content per unit mass of collagen, and ρ c,tot signifies the total collagen content as a proportion of the tissue's total wet weight (WW).
To ascertain water content, the initial step involves measuring the tissue sample's wet weight (WW) (Huyghe et al., 2003;Malandrino et al., 2015;Ruiz et al., 2016), followed by lyophilization to obtain the dry weight (DW).These measurements facilitate the calculation of the initial total water content (n f,0 ) and, subsequently, the initial solid fraction and the current fluid fraction: Thus, using n f,0 as a foundational value, the equations seamlessly connect n s,0 and n f with the IVD deformations represented by J, establishing a coherent framework for relating IVD composition measurements to mechanical modeling parameters.Then, the total fluid volume ratio is calculated using the void volume in the medium (dV v ) and the total volume of the medium (dV ): To estimate the proteoglycan and total collagen contents, previous works propose digesting the dried samples in a papain solution.The digested solutions were then used (i) to determine the content of sulfated glycosaminoglycans (sGAG) through a dimethyl methylene blue (DMMB) assay (Farndale et al., 1986) and (ii) to achieve a measure for collagen content according to hydroxyproline measurements through the chloramine-T assay (Huszar et al., 1980).
We calculated the initial fixed charge density (c f,0 ) per total hydrated tissue volume, from which c f (Equation S16) is derived to be dependent of J, using the expression (Narmoneva et al., 1999): Here, z cs , MW cs , and c cs are the valency (2 mEq/mmol), the molecular weight (513,000 µg/mmol), and the concentration (in µg/mL) of chondroitin sulfate, respectively.The sGAG content measured through the DMMB assay is assumed to be equivalent to the chondroitin sulfate content, i.e., c cs is the amount of sGAG divided by the sample's water content.To obtain ρ c,tot , in Equation S24, the initial collagen content (µg/mg DW) was estimated from hydroxyproline content by using 7.6 as the mass ratio of collagen to hydroxyproline (Sivan et al., 2006): See the works of Huyghe et al. (2003); Malandrino et al. (2015); Ruiz et al. (2013Ruiz et al. ( , 2016Ruiz et al. ( , 2018) ) for more information on the evolution of the used model.

Permeability
The tissue's AF and NP hydraulic permeability (κ) are strain-dependent according to the following expression as Wilson et al. (2006a) developed and then adapted by (Schroeder et al., 2007) for the IVD: where α stands for the initial permeability at zero strain, and M is a positive constant that governs volumetric strain dependency, and n f,exf is calculated in Equation S17.
The CEP is also strain-dependent, following this equation (Mow et al., 1980;Argoubi and Shirazi-Adl, 1996): where e is the void ratio, which is the ratio of the current pore volume (i.e., fluid) to the current volume of the solid matrix, and e 0 is the initial void ratio.The void ratio is related to the initial and current water content/porosity of the tissue, n f,0 and n f , according to the following expression:
S27) Figures corresponding to each of the volumes.Posterior Transition Zone (PTZ), CNP: Center of the Nucleus, Pulposus, and ATZ: Anterior Transition Zone, and surfaces CT: Cartilage Endplate Top, CB: Cartilage Endplate Bottom, used to evaluate the morphological factors on the mechanical variables.

Supplementary Material: Unveiling IVD's Interaction Between Morphology and Mechanics CNP Step
Figures corresponding to each of the volumes, Posterior Transition Zone (PTZ), CNP: Center of the Nucleus, Pulposus, and ATZ: Anterior Transition Zone, and surfaces CT: Cartilage Endplate Top, CB: Cartilage Endplate Bottom, used to evaluate how big is the impact of the morphological factors on the mechanical variables.